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ABSTRACT 

Using a collapsar progenitor model of MacFadyen & Woosley we have simulated the propagation of an ax- 
isymmetric jet through a collapsing rotating massive star with the GENESIS multi-dimensional relativistic hy- 
drodynamic code. The jet forms as a consequence of an assumed (constant or variable) energy deposition in the 
range 10^" erg s"' to 10^' erg s"' within a 30° cone around the rotation axis. The jet flow is strongly beamed 
few degrees), spatially inhomogeneous, and time dependent. The jet reaches the surface of the stellar progenitor 
= 2.98 X 10'"cm) intact. At breakout the maximum Lorentz factor of the jet flow is 33. After breakout the 
jet accelerates into the circumstellar medium, whose density is assumed to decrease exponentially and then being 
constant pext = 10"^ gcm"^. Outside the star the flow begins to expand also laterally (v ^ c), but the beam remains 
very well collimated. At a distance of 2.54R^, where the simulation ends, the Lorentz factor has increased to 44. 

Subject headings: Hydrodynamics — Methods: Numerical — Relativity — Gamma Rays: Bursts, Theory 



1. INTRODUCTION 

Catastrophic collapse events have been proposed to explain 
the energies released in a gamma-ray burst (GRB) includ- 
ing mergers of compact binaries (Pacynski 1986; Goodman 
1986; Eichler et aJ. 1989; Mochkovitch et aJ. 1993), collapsars 
(Woosley 1993) and hypernovae (Pacynski 1998). According 
to the current view these models require a stellar mass black 
hole (BH) which accretes up to several solar masses of mat- 
ter powering a pair fireball. If the baryon load of the fire- 
ball is not too large, baryons are accelerated together with 
e''"e" pairs to Lorentz factors > 10^ (Cavallo & Rees 1978). 
Such relativistic flows are supported by radio observations of 
GRB 980425 (Kulkarni et al. 1998b). Sphericafly symmetric 
fireballs have been studied by several authors by means of 
ID Lagrangian hydrodynamic simulations (e.g., Panaitescu et 
al. 1997; Panaitescu & Meszaros 1998; Kobayashi, Piran & 
Sari 1999). Recently, it has been argued that the rapid tem- 
poral decay of several GRB afterglows is more consistent with 
the evolution of a relativistic jet after it slows down and spreads 
laterally than with a spherical blast wave (Sari, Piran & Halpern 
1999; Halpern et al. 1999; Kulkarni et al. 1999; Rhoads 1999). 
The lack of a radio afterglow in GRB 990123 provides indepen- 
dent evidence for jet-like geometry (Kulkarni et al. 1999b). 

2. INITIAL MODEL AND NUMERICAL SETUP 



MacFadyen & Woosley (1999; MW99) have explored the 
evolution of rotating helium stars (M^^IOM©) whose iron 
core collapse does not produce a successful outgoing shock, but 
instead forms a BH surrounded by a compact accretion torus. 
Assuming an enhanced efficiency of energy deposition in polar 
regions MW99 obtain relativistic jets along the rotation axis, 
which are highly focused and seem to be capable of penetrat- 
ing the star However, as their simulations are Newtonian, they 
obtain flow speeds which are superluminal. 

We have performed axisymmetric relativistic simulations us- 
ing a 14M0 collapsar model from MW99. When the central 
BH has acquired a mass of 3.162Mq we map the model to 
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our computational grid. In a consistent collapsar model a jet 
will be launched by any process which gives rise to a local de- 
position of energy and/or momentum, as e.g., i/j>-annihilation, 
or magneto-hydrodynamic processes. We mimic such a pro- 
cess by depositing energy at a prescribed rate homogeneously 
within a 30° cone around the rotation axis. In the radial direc- 
tion the deposition region extends from the inner grid bound- 
ary located at 200 km to a radius of 600km. We have inves- 
tigated constant energy deposition rates E = 10^" ergs"', and 
E = 10^' erg s"', and a varying deposition rate with a mean value 
of 10^" erg s"' . The constant rates roughly bracket the expected 
E of collapsar models, while the varying rate mimics, e.g., time- 
dependent mass accretion rates resulting in time-dependent vv- 
annihilation (MW99). 

The simulations were performed with the multidimensional 
relativistic hydrodynamic code GENESIS (Aloy etaJ. 1999) us- 
ing a 2D spherical grid with 200 radial zones spaced logarithmi- 
cally between the inner boundary and the surface of the helium 
star at R^, = 2.98 x 10"^ cm. Assuming equatorial symmetry we 
performed simulations with 2°, 1° and 0.5° angular resolution. 
In the latter case the grid consists of 60 uniform zones covering 
the polar region (0° < 6 < 30°) and 40 nonuniform zones log- 
arithmically distributed between 30° < 9 < 90°. The gravita- 
tional field of the BH is described by the Schwarzschild metric. 
Effects due to the self-gravity of the star are neglected, i.e., we 
consider only the gravitational potential of the BH. The equa- 
tion of state includes non-relativistic nucleons treated as a mix- 
ture of Boltzmann gases, radiation, and an approximate correc- 
tion due to e^e"-pairs as described in Witti, Janka & Takahashi 
(1994). Complete ionization is assumed, and the effects due to 
degeneracy are neglected. We advect nine non-reacting nuclear 
species which are present in the initial model: C'^, O'^, Ne^", 
Mg^"*, Si^^, Ni^*, He"*, neutrons and protons. 



3. RESULTS 

3.1. Constant small energy deposition rate (Model C50) 
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For a constant E = 10^*' erg s~' a relativistic jet forms within 
a fraction of a second and starts to propagate along the rota- 
tion axis with a mean speed of 7.8 x 10^ cms" ^ (Fig. 1). The 
jet exhibits all the morphological elements of the Blandford & 
Rees (1974) jet model in the context of classical double radio 
sources: a terminal bow shock, a narrow cocoon, a contact dis- 
continuity separating stellar and jet matter, and a hot spot. Fig. 1 
shows that the density structure of the star does not change no- 
ticeably during the whole evolution. This, a posteriori, justifies 
our neglect of the self-gravity of the star 

The propagation of the jet is unsteady, because of density 
inhomogeneities in the star. The Lorentz factor of the jet, F, in- 
creases non-monotonically with time, while the density drops 
to ~ 10"*gcm"^ (Fig. 2). The density profile shows large vari- 
ations (up to a factor of 100) due to internal shock waves. The 
mean density in the jet is ~ 10"^ gcm"^. Some of the internal 
biconical shocks, which develop during the jet's propagation, 
recolUmate the beam. They may provide the "internal shocks" 
proposed to explain the observed gamma-ray emission. A par- 
ticularly strong recollimation shock wave (hardly evident at low 
resolution) forms early in the evolution. A very strong rarefac- 
tion wave behind this recolUmation shock causes the largest lo- 
cal acceleration of the beam material giving rise to a maximum 
in the Lorentz factor. When the jet encounters a region along 
the axis where the density gradient is positive (at log r « 8 . 1 and 
log r f» 8.6) the jet's head is decelerated, while a central channel 
in the beam is cleaned by outflow into the cocoon through the 
head, which accelerates the beam. The combination of both ef- 
fects (deceleration of the head and beam acceleration) increases 
the strength of the internal shocks. Within the jet the mean 
value of the specific internal energy, ~ 10^"^- 10^^ ergg"\ or 
~ O(c^). The mean temperature ~ 5 x 10^ K (well below the 
pair creation threshold) implying that the pressure is radiation 
dominated in accordance with our simphfied EOS. 

The relativistic treatment of the hydrodynamics leads to a 
quaUtatively similar (formation of a jet), but quantitatively very 
different evolution than in MW99. According to their Fig. 27 
the jet propagates 7 000 km within the first 0.82 s. Furthermore, 
MW99 infer an asymptotic F ~ 10, and find a half opening an- 
gle, f2, for their jet of ~ 10°. In our simulation, at the same time 
for the same angular resolution 2°) and E the head reaches 
a radius of 30000 km, but the maximum Lorentz factor (Fmax) 
is only 4.62 at ^ 12200 km. Such quantitative difference is ex- 
pected, in part, due to different mapping time and inner bound- 
ary radius between the two calculations. Initially, in our simu- 
lations Vl is between 6° to 8° depending on angular resolution. 
At ~ 1.5s the strong recollimation shock reduces fi^l". 

We find that some results strongly depend on angular reso- 
lution, the minimum acceptable one being 0.5° (at least near 
the axis). The morphology of the jet is richer at higher reso- 
lution. At 0.5° angular resolution Tmux ~ 15-20 at a radius 
~ 8 X 10^ cm at jet breakout. Within the uncertainties of the jet 
mass determination due to finite zoning and the lack of a pre- 
cise numerical criterion to identify jet matter, the baryon load 
{rj = Mc^/Edg^os with fidepos = / Edt) decreases with increas- 
ing resolution. In the highest resolution run we find an average 
baryon load of 77 ~ 1.3 at jet breakout (see also Sect. 4). 

3.2. Constant large energy deposition rate (Model C51 ) 

Enhancing £ by a factor of ten (to 10^' erg s"' ), the jet flow 
reaches larger Lorentz factors. We observe transients during 
which the Lorentz factor becomes as large as 40. After 1 .2 s the 



Lorentz factor steadily increases from 22 to 33. The jet propa- 
gates faster than in model C50. The time required to reach the 
surface of the star is 2.27 s instead of 3.35 s. At breakout the jet 
is less colUmated (O ~ 10°). The strong recollimation shock 
present in model C50 is not so evident here. Instead, several 
biconical shocks are observed within a very knotty beam and 
the Lorentz factor near the head of the jet is larger (~ 22 in the 
final model) because, due to the larger E, the central funnel is 
evacuated faster, and because the mean density of the jet is 5 
times smaller than in model C50 (ry being twice as large). 

3.3. Varying energy deposition rate (Model V50) 

We have computed a model where the mean energy depo- 
sition rate (lO^'^ergs"') randomly varies on time scales of a 
few milUseconds and the amplitude by a factor of ten. Com- 
pared to model C50 the jet structure is more knotty and also 
richer in shocks, particularly inside the first 10^ cm (where the 
radial resolution is large enough to capture the finest structures 
imprinted on the flow by the time variability of the deposition 
rate). At breakout Fmax = 26.81, which is almost twice as large 
as the one found in model C50. Thus, a variable E is more 
efficient in converting internal energy into kinetic energy, and 
in this case, the internal shocks are stronger and more numer- 
ous. The mean propagation speed is similar in both models, 
although the instantaneous velocity of the jet's head is clearly 
different. Behind the strongest recoUimation shock f2 < 1° in 
both models. 

3.4. Evolution after jet breakout 

The structure of the circumsteUar medium will influence the 
characteristics of the GRB and of the subsequent afterglow. 
Thus, a continuation of the simulations beyond jet breakout 
is necessary. In order to satisfy the conditions for accelerat- 
ing shocks (Shapiro 1979) we endowed the star with a Gaus- 
sian atmosphere, which at /?a = 1-8/?* passes over into an ex- 
ternal uniform medium with a density 10"^ gcm"^ and a pres- 
sure 10~**/?(/?*). The computational domain is extended to 
Rt = 2.54/J, with 70 additional zones. The evolution after jet 
breakout has been computed for models C50 and C51. In both 
cases the jet reaches Ri after ^ 1 .8 s (measured from breakout). 
Its mean propagation velocity is ~ 0.85 c, which is almost three 
times faster than the velocity of the head inside the star (0.30 c 
in model C50; 0.44 c in model C5 1). 

The evolution after jet breakout consists of three distinct 
epochs (Fig. 3). The first one lasting 0.35 s is characterized by 
a head velocity of 0.48 c and a small sideways expansion. Dur- 
ing the second phase (of 0.3 s) the jet head accelerates to 0.91 c, 
because of the steep external density gradient, and because the 
flux of axial momentum is still important compared to pressure. 
The sideways expansion is still sub-relativistic (« 0.008 c), and 
f2 of the beam increases to « 10°. During the final 1.2 s the 
bow-shock propagates within the uniform part of the ambient 
medium leading to a rapid (P ~ 5) lateral spreading (Fig. 3). 

The shape of the expanding bubble is prolate (Figs. 1 and 
3) during the post-breakout evolution. However, when the jet 
reaches the uniform part of the circumsteUar environment, the 
bubble widens due to the faster sideways expansion. We expect 
a more isotropic expansion when most of the bubble is inside 
the uniform medium, and when it is pressure driven (in particu- 
lar if the energy deposition is switched off). The Lorentz factor 
near the boundary of the cavity blown by the jet grows from 1 
(at jet breakout) to ~ 3 in both models decreasing with latitude. 
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At the end of the simulation r^ax is 29.35 (44.17) for model 
C50 (C51), which is still smaller than the ones required for the 
fireball model (Piran 1999). 

4. DISCUSSION AND CONCLUSIONS 

Energy deposition in the polar regions of a collapsar model 
gives rise to both the formation and propagation of a relativistic 
jet through the mantle and envelope of the star, and to a super- 
nova explosion. The jet has a small opening angle 8°) and 
possesses a highly collimated 1°), ultra-relativistic core in 
which the Lorentz factor reaches a value of r^ax = 44 (model 
C5 1 ) at the end of the simulation about 2 s after shock break- 
out. The equivalent isotropic kinetic energy (see MacFadyen, 
Woosley & Heger 1999) slightly exceeds 10^^ erg for model 
C51 (10^^ erg for model C50) within2° (5°) of the rotation axis 
dropping by a factor of 10 within 17° (10°). The inner region 
contains 8 x 10~'^Mq with < T >^ 4 (Fig. 1). For a larger E, 
the jet and in particular the cocoon are less collimated, because 
a harder driven jet also expands stronger laterally. 

The rest-mass density and the internal energy strongly vary 
in space and time within the jet giving rise to a very inhomo- 
geneous baryon load, i.e., the concept of 77 as a global param- 
eter is useless. Instead it is more appropriate to discuss the 
efficiency of energy conversion in terms of incremental baryon 
loads considering only matter within a given range of F-values. 
Although we find an average baryon load of the jet of 77 ~ 1, 
some parts of the flow have a baryon load as low as ~ 10"^ or 
even less. After jet breakout fj decreases by a factor 4 in less 
than 1.8 s. If this trend continues even fj ^ lO""* within 9 s. 

In model C51, at the end of the simulation, ~ 2Mq have 
a Lorentz factor of less than three, 3 x 1Q~^Mq move with 
3 < r < 10, and for 2 x 10"^ flie Lorentz factor L > 10 
(Fig. 4). The latter two masses reduce to 2 x 10~^ Mq and 
2 X 10"^ M0 for model C50. Except for the very early evolution 
(f < 1 s) the amount of matter moving at moderate (3 < T < 10) 
and highly (F > 10) relativistic velocities increases by a factor 
~ 3 every second, i.e., if the central engine is active for another 
5 s, at the assumed energy deposition rate, ^ IQ~* Mq will move 
with r > 10. As Fmax is also rapidly increasing, Lorentz factors 
of several hundreds might be reached before the central engine 
is switched off. For models where the released total energy is 



equal (C50 and V50) Fmax is higher (by a factor of two) for a 
time-dependent E. Determining the efficiency of energy con- 
version £ = Ek/Eispos is hampered by the fact, that the kinetic 
energy Eit of a relativistic fluid is not a well defined quantity. If 
Ek = J pF(F- l)dV, we find £ = 1 .6 for model C50 (and f = 2.8 
for model C51) at the end of the simulation. These efficien- 
cies are obtained considering only matter with radial velocities 
> 0.3c and specific internal energy densities > 5 x 10''' ergg"', 
i.e., matter in the jet. Note that efficiencies larger than one can 
arise, because (i) there are also other large sources of energy 
(e.g., gravitational, internal) available, and (ii) because matter 
is entrained into the jet which does not originate from the de- 
position region. Thus, efficiencies larger than one suggest that 
the local energy deposition is efficient triggering conversion of 
other forms of energy into kinetic energy. 

In our simulations the jet reaches the stellar surface intact 
(propagating over three decades in radius). This result may also 
hold for other less specific initial conditions. A more spherical 
density stratification might decrease the collimation of the jet, 
but the outflow might also be initiated mostly by momentum 
deposition instead of pure energy deposition (e.g., by MHD ef- 
fects). The propagation of the jet after breakout will depend on 
the density stratification of the circumstellar medium. Thus, 
further simulations with different environments are planned. 
We note in this respect that the post-breakout propagation is 
similar in models C50 and C5 1 suggesting that a lower value of 
Pext will not change the dynamics quantitatively. In our models 
the jet has only reached a radius of 7.5 x 10'° cm at the end of 
our simulations, which is 10^ to 10"* times smaller than the dis- 
tance at which the fireball becomes optically thin. Determining 
whether a GRB will eventually be produced requires to com- 
pute the further evolution of the jet. As the jet has stayed col- 
limated in the star it might remain focused over the next three 
decades, too. 

This work has been supported in part by the Spanish DGES 
(grant PB97-1432) and the CSIC. MAA expresses his gratitude 
to the Conselleria d'Educacio i Ciencia de la Generalitat Valen- 
ciana for a fellowship. We would like to thank Stan Woosley 
for his enthusiasm in promoting this work, and Thomas Janka 
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Fig. 1. — Contour maps of the logarithm of the rest-mass density after 3.87 s and 5.24 s (left two panels), and of the Lorentz factor 
(right panel) after 5.24 s. X and Y axis measure distance in centimeters. Dashed and solid arcs mark the stellar surface and the 
outer edge of the exponential atmosphere, respectively. The other solid line encloses matter whose radial velocity > 0.3c, and whose 
specific internal energy density > 5 x lO'^ erg g~' . 
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Fig. 2. — Rest-mass density (top) and Lorentz factor (bottom) vs radius along the symmetry axis for model C50 at ? = Os (long 
dashed), ? = 1 .44 s (dotted), ? = 3.87 s (dashed) and ? = 5.24 s (solid), respectively. 
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Fig. 3.— Axial and lateral size of the jet driven bubble versus time since jet breakout for models C50 and C5 1 . 
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Fig. 4.— Mass (binned in AF = 0.5 intervals) versus Lorentz factor for model C50 at shock breakout (dash-dotted) and at ? = 5.24 s 
(sohd), respectively. 



